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ABSTRACT 

Magnetic reconnection is thought to be the driver for many explosive phenomena in the universe. 

The energy release and particle acceleration during reconnection have been proposed as a mechanism 
for producing high-energy emissions and cosmic rays. We carry out two- and three-dimensional kinetic 
simulations to investigate relativistic magnetic reconnection and the associated particle acceleration. 

The simulations focus on electron-positron plasmas starting with a magnetically dominated, force-free 
current sheet (a = /{AirnerrieC^) ^ 1). For this limit, we demonstrate that relativistic reconnection 

is highly efficient at accelerating particles through a first-order Fermi process accomplished by the 
curvature drift of particles along the electric field induced by the relativistic flows. This mechanism 
gives rise to the formation of hard power-law spectra / oc (7 — 1 )“^ and approaches p = 1 for 
sufficiently large a and system size. Eventually most of the available magnetic free energy is converted 
into nonthermal particle kinetic energy. An analytic model is presented to explain the key results 
and predict a general condition for the formation of power-law distributions. The development of 
reconnection in these regimes leads to relativistic inflow and outflow speeds and enhanced reconnection 
rates relative to non-relativistic regimes. In the three-dimensional simulation, the interplay between 
secondary kink and tearing instabilities leads to strong magnetic turbulence, but does not significantly 
change the energy conversion, reconnection rate, or particle acceleration. This study suggests that 
relativistic reconnection sites are strong sources of nonthermal particles, which may have important 
implications to a variety of high-energy astrophysical problems. 

Subject headings: acceleration of particles - magnetic reconnection - relativistic process - gamma-ray 
bursts: general - galaxies: jets - pulsars: general 


1 . INTRODUCTION 

Magnetic reconnection is a fundamental plasma pro¬ 
cess that rapidly rearranges magnetic topology and con¬ 
verts magnetic energy into various forms of plasma ki¬ 
netic energy, including bulk plas ma flow, thermal and 


nonthermal plasma distributions (Kulsrud 1998 Priest 


fc Forbes||2QQQ ). It is thought to play an irnportant role 


during explosive energy release processes of a wide va¬ 
riety of laboratory, space and astrophysical systems in¬ 
cluding tokamak, planetary magnetospheres, solar flares, 
and high-energy astrophysical objects. Relativistic mag¬ 
netic reconnection is often invoked to explain high-energy 
emissions and ultra-high-energy cosmic ray s from ob 


Arons||2012 Uzdensky & Spitkovskyi2014|), jets from ac- 

tive galactic nuclei (AGJN 

; |de Gouveia dal Pino & Lazar- 

ian 2005 

IGiannios et al. 

2 OO 9 I), and gamma-ray bursts 

( GRHs; iThompson 1994 

Zhang Van 

2011 McKinney 

& Uzdensky||2012p. In those systems, t 
parameter cr = /{AirnemeC^)^ is oftei 

much larger than unity a ^ 1 and 1 
approaches the speed of light va ^ c. 
observed high-energy emissions, often a 
conversion mechanism is required (e.g.. 

he magnetization 

1 estimated to be 
ffie Alfven speed 
To explain the 
n efficient energy 

Zhang et al.|2007 

Celotti & Ghisellini 2008 

Zhang & Yan 

2011t|Zhang et al. 
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plasma flow energy into thermal and nonthermal ener¬ 
gies in low-cr flows, are inefficient in dissipating magneti¬ 
cally dominated flows, where most of the energy is stored 
in magnetic fields. In these regimes, magnetic reconnec¬ 
tion is the primary candidate for dissipating and convert¬ 
ing magnetic energy into relativistic particles and subse¬ 
quent radiation. Understanding magnetic reconn ection 
is also important to solve the so-called a-problem dCoro- 
niti||199Q| ILyubarsky fc Kirk||2QQf| [Kirk fc !Skjaeraas~ 
2QQ3[ Forth et al.|2Q13p , where strong magnetic dissipa- 


tion may be required to convert the magnetically domi¬ 
nated flow (cr 1 ) to a matter-dominated flow (a ^ 1 ). 
However, the detailed physics of relativistic magnetic re¬ 
connection, including the magnetic reconnection rate, en¬ 
ergy conversion and particle acceleration, are not well 
understood. 


¥ 

the 


1994) and Lyutikov & Uzdensky 


jBlac kman & Field 

(|2UU3|) have studied the properties of relativistic mag- 
netic reconnection using the extended Sweet-Parker and 
Petschek models. They found that when a ^ 1 the out¬ 
flow speed Uout approaches the speed of light, and the 
rate of relativistic magnetic reconnection and inflow ve¬ 
locity Uin may increase compared to the nonrelativistic 
case. This is because of the enhanced outflow density 
arising from the Lorentz contraction of plasma passing 
through the diffusion region Uin oc Uout^out/^in^ where 
Tout and F^^ are Lorentz factors of outflo ws and in¬ 
flows , respectively. However, later analysis (Lyubarsky 


2005) showed that for a pressure-balanced current layer 
(H’^/Stt ^ nk{Ti + Tg)), the thermal pressure constrains 
the outflow speed to be mildly relativistic and hence the 
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effect of Lorentz contraction is negligible. Although the 
rate of relativistic magnetic reconnection is reported to 
enhance in a number of studies using different numeri- 
cal m odels dZenhani et al.|2QQ9 Bessho & Bhattacharjee 


2014). For a Harris current layer, it was found that a 


20121 I'Fakamotol 120131 Ibironi fc bpitkovsl^ |2014| |Guo 

et al.jpQMf Melzani et al.||2014ap , its nature is not clear. 

4'his issue is recently revisited by carefully analyzing re- 
sults from fully k inetic two-dimensional (2D) simulations 
(Liu et al. 2015), which shows that the plasma density 
and pressure around the X-line drop significantly as the 
initial high-pressure region is depleted during reconnec¬ 
tion. This results in a reconnection region with a ^ 1 
and a relativistic inflow speed Vin ^ c. The local recon¬ 
nection rate across the diffusion region is well predicted 
by a simple model that includes the Lorentz contrac¬ 
tion. However, the extension of these results to three- 
dimensional (3D) kinetic simulations was not considered. 

Plasma energization during magnetic reconnection has 
been extensively discussed in the literature. However, 


the primary acceleration m echanism is still unclear. Ro¬ 


manova & Lovelace (1992) analyzed particle motions in a 
large-scale reconnection region and predicted a spectrum 
dN/dj = 7 ~^ w ith p = 1.5 for the pair plasma case. 


Litvinenko (1999) found a solution with a spectral index 
p = 2 when particles are accelerated in a direct elec¬ 
tric field associated with magnetic reconnection. Using a 
model for the motio ns of particles in a stea dy magnetic 
reconnection region, Larrabee et al. (2003) have found 
strong particle acceleration in the reconnection layer and 
obtained a hard energy spectrum with a spectral index 
of about p = 1. The first-order Fermi acceleration in 
converging reconnection inflows has been discussed ( de 
Gouveia dal Pino fc: Lazarian||2005} I Lazarian fc: Ophe r 


200^I^wal et al.|2012p . Drury 


studied the accel- 

eration in a reconnection layer including energy change 
in both inflows and outflows and show fluid compression 
is crucial for efficient particle acceleration. Test-particle 
approach has been applied to interpret the strong parti¬ 
cle acce leration responsible for y-ray flares from the Grab 


pulsar ( Gerutti et al.||2Q12 ). The results show that mag¬ 


netic reconnection may be the site of extreme particle 
acceleration required to expl ain high energy em issions 
from the Grab flares (see also Gerutti et al. 2013). Self- 


consistent kinetic simulations have been widely used to 
study plasma dynamics and particle energization during 
magnetic reconnection. Most of previous kinetic studies 
have focused on the regime with cr < 1 , and found a num¬ 
ber of acceleration mechanisms such as direct accelera- 


tion at X-line regions ([Drake et al.[[2005| Fu et al. 

[2006 

Pritchett[[2006 Oka et al. 2010| [Huang et al.|[201U 

I and 

fermi-type accel 
flows within maj 

eration in reconnection induced p 

lasma 

genetic islands ([Drake et al. 20061 

2010 

Oka et al.[[2010 Huang et al.[[ 2010 [). The high-a regime 


the Harris equil i brium (Zenitani fc Hoshino| 2QQ1[ 2007 
Liu et al. 201 1| Bessho fc Bhattacharjee|| 2012 [ |(Jerum 
et al.| ^13[ ISironi fc Spitkovsky 2014| Melzani et al. 

2014b i Werner et al. 2014p . However, the initial con- 


dition employed in these studies requires a hot plasma 
component inside the current sheet to maintain force 
balance, which may not be justified for high-cr plasmas. 
Recently, several studies have repo rted hard power-law 
distributions 1 < p < 2 when a ^ 1 dSironi fc Spitkovsky 
|2014[[Guo et al.|2014 Melzani et al.|2014b[| Werner et a l. 


power-law distribution can be obtained by subtracting 
the initial hot plasnia component in the current layer 


(ISironi fc S pitkovsky||201 4| [Melzani et al.|j2 014b| [Werner 

et al.||2014|). In contrast, [Guo et al.|(|2014|) used a torce- 


et al.|[2014). In contrast, Guo et al. (|20l4|) used a force- 

tree current sheet that does not require tne hot plasma 
population and showed that the energy distribution of 
particles within the entire reconnection layer develops a 
power-law distribution. In this study, the primary accel¬ 
eration mechanism was demonstrated to be a first-order 
Fermi mechanism resulting from the curvature drift of 
particles in the direction of the electric field induced by 
the relativistic flows. This mechanism gives rise to the 
formation of hard power-law spectra / oc (7 — 1 )“^ with 
spectral index approaching p = 1 for a sufficiently high 
a and a large system size. An analytical model was de¬ 
veloped to describe the main feature of the simulations 
and it gives a general condition for the formation of the 
power-law particle energy distribution. The solution also 
appears to explain simulations from the Harris current 
layer, in which the particles initially in the current layer 
form a heated thermal distribution and particles injected 
from the upstre am region are accelerated into a power- 
law distribution dSironi fc S pitkovsky 2014 Melzani et al. 
2014b| I Werner et al.||2U14|). ' 

Another important issue is the influence of 3D dynam¬ 
ics that may significantly modify the reconnection rate, 
energy release, and particle acceleration process. Re¬ 
cently, the rate of 3D nonrelativistic magnetic reconnec¬ 
tion has been explored and compared with 2 D simula¬ 
tions in a number of non-re lativistic studies (Liu et al. 


2013 Daughton et al. 2014), which showed only mod¬ 


est difterences between 2 D and 3D simulations although 
strong 3D effects emerge as the tearing mode dev elops 
over a range of oblique angles (Daughton et al. 2011). For 
relativistic magnetic r econn ection with a pair plasma, 
Sironi & Spitkovsky (2014) reported a factor of four 
decrease of reconnection rate for 3D simulations com¬ 
pared to 2 D simulations. This is in contrast to [Guo 
et al. (2014), who observed similar reconnection rate and 


energy conversion betw een 2 D and 3D simulations, al¬ 
though the kink mode (Daughton 1999) strongly inter¬ 
acts with the tearing mode lea ding to a turbulent recon¬ 
nection layer (Yin et al.[[2008). For particle acceleration 
in 3D reconnection simulations, early studies reported 
that the drift kink instability can modify the electric and 
magnetic field structures in an antiparallel r econnection 
layer and prohibit nonther mal acceleration (Zenitani & 
Hoshino[[2005b| [2007| [2008[). However, recent large-scale 
3D simulations have found strong nonthermal particle 
spectra are produced even when the kink mode is active 


iLiu et ah 2011 Sironi & Spitkovsky 2014 Guo et al 


014[ ). Earlier large-scale 3D studies have shown th e de^ 


velopment of turbulence in the reconnection layer (Yin 


et al. 


2008), but its effect to energetic particle acceler 


ation IS unknown. It is therefore important to further 
study particle acceleration in relativistic regimes using 
large-scale 3D kinetic simulations. 

In this paper, we perform 2D and 3D fully kinetic simu¬ 
lations starting from a force-free current sheet with uni¬ 
form plasma density and temperature to model recon¬ 
nection over a broad range in the magnetization param¬ 
eter a — 0.25 - 1600. This paper builds upon earlier 
work ( Guo et ^[2014 ) and gives further details regard- 
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ing the plasma dynamics and particle acceleration during 
relativistic magnetic reconnection in the high-a regime. 
We also present detailed results from a 3D simulation 
that shows a turbulent reconnection layer arising from 
the interaction between the secondary tearing and kink 
modes. In Section 2 , we describe the numerical methods 
and parameters. Section 3 discusses the main results of 
the paper. In section 4, we present an analytical model 
that explains the main feature of particle acceleration in 
the simulations. The implications from this work for a 
range of astrophysical problems are discussed in Section 
5 and our conclusions are summarized in Section 6 . In 
addition, we have also explicitly examined the numerical 
convergence for this problem and the effect of numeri¬ 
cal heating in our simulations, which is discussed in the 
Appendix. 


2 . NUMERICAL METHODS 

We envision a situation where intense current sheets 
are developed within a magnetically dominated plasma. 
Earlier work in non-relativistic \ow-/3 plasmas has shown 
that the gradual evolution of the magnetic field can lead 
to formation of intense nearly force-free current layers 
where magnetic reconnection may b e triggered (Titov 
et al.| |2003t |Galsgaard et al.j |2003|) . In the present 
study, the critical parameter s the magnetization pa¬ 
rameter defined as a = /{AirnerrieC^)^ which roughly 
corresponds to the available magnetic energy per par¬ 
ticle. The numerical simulations presented in this pa¬ 
per are initialized from a force-free c urrent layer with 
B = B()tanh(z/X)x ^ B()sech(z/X)y ( |Che et 2011 


Liu et al. 2013, 2014), which corresponds to a mag- 


netic held with magnitude Bq rotating by 180° across 
the central layer with a half-thickness of A. No exter¬ 
nal guide field is included in this study but there is 
an intrinsic guide field By associated with the central 
sheet. The plasma consists of electron-positron pairs 
with mass ratio mijm^ = 1. The initial distributions 
are Maxwellian with a spatially uniform density no and 
a thermal temperature {kTi = kTe = 0.36meC^). Par¬ 
ticles in the central sheet have a net drift = — Ue 
to represent a current density J = eno{\Ji — Ug) that 
is consistent with V x B = 47 rJ/c. Since the force-free 
current sheet does not require a hot plasma component 
to balance the Lorentz force, this initial setup is more 
suitable to study reconnection in low (3 and/or high-cr 
plasmas. The full p article simulations a re performed us- 
i ng the VPIC code (iBowers et al.||2QQ9|) and NPIC cod e 
(jDaughton et al.||2006t |D aught on Karimabadi||2UU7|), 
both of which solve Maxwell equations and push particles 
using relativistic approaches. The VPIC code directly 
solves electric and magnetic fields in Maxwell equations, 
whereas in the NPIC code the fields are advanced us¬ 
ing the scalar and vector potentials. Therefore the two 
codes have different algorithms and can be used to cross¬ 
check numerical results for this problem. We find that 
they give consistent results for this problem. In addition, 
we have developed a particle-tracking module to analyze 
the detailed physics of the particle energization process. 
In the simulations, we define and adjust a by changing 
the ratio of the electron gyrofrequency Ugg = eBjim^c) 
to the electron plasma frequency cjpe = jme^ 

(7 = B ^/= (Uge/co’pe)^- For 2D simulations. 


we have performed simulations with a = 0.25 ^ 1600 
and box sizes 'x Lz = SOOdi x 194(i^, GOOdi x 388di^ 
and 1200di x 776d^, where di is the inertial length cjujye. 
Eor 3D simulations, the largest case is x Ly x Lz = 
SOOdi X 194:di X SOOdi with a = 100. Eor high-cr cases 
(cr > 25), we choose cell sizes Ax = Ay = 1.46/v^d^ 
and Az = 0.95/ ^/o^di, so the particle gyromotion scale 
^ "^the/V^di is resolved. The time step is chosen to 
correspond to a Courant number Cr = cAt/Ar = 0.7, 
where Ar = Ax Ay Az/{Ax Ay -f AyAz + AxAz). The 
half-thickness of the current sheet is A = Gdi for cr < 100, 
12di for cr = 400, and 2Adi for cr = 1600 in order to sat¬ 
isfy the drift velocity < c. Eor both 2D and 3D sim¬ 
ulations, we have more than 100 electron-positron pairs 
in each cell. The boundary conditions for 2D simula¬ 
tions are periodic for both fields and particles in the 
x-direction, while in the ^-direction the boundaries are 
conducting for the field and reflecting for the particles. 
In the 3D simulations, the boundary conditions are peri¬ 
odic for both fields and particles in the ^-direction, while 
the boundary conditions in the x and z directions are the 
same a s the 2 D cases. A weak long-wavelength pertur¬ 
bation (|Birn et al.| 2001 |) with Bz = O.OSBq is included to 
initiate reconnection. I'he parameters for different runs 
are summarized in Table 1, which also lists key results 
such as maximum energy of particles, spectral index, the 
fraction of kinetic energy converted from the magnetic 
energy and the portion of energy gain arising from the 
perpendicular electric fields. 

Using the set of numerical parameters described above, 
all of the simulations show excellent energy conservation 
with violation of energy conservation less than 10“^ of 
the total energy in all cases. However, we note that to 
accurately determine the particle energy spectra, the vi¬ 
olation in energy conservation should be smaller than the 
initial plasma kinetic energy, which is only a small frac¬ 
tion of the total energy for the problem we study. Cau¬ 
tion is needed when using a small number of particles 
per cell and a small initial plasma kine tic energy in the 
simulations (|Sironi & Spitkovsky||2014|), since numerical 
heating may significantly modify the particle distribu¬ 
tion. In the Appendix, we have extensively tested how 
the numerical convergence varies with the initial plasma 
temperature, cell size, number of particles per cell, and 
time step. Eor all the cases we present in the main paper, 
the violation of energy conservation is a few percent of 
the initial kinetic energy in the system, meaning effects 
such as numerical heating have a negligible influence on 
the simulated energy spectra. 


3. SIMULATION RESULTS 
3.1. General feature and energy eonversion 

Eigure 1 gives an overview of the evolution of the cur¬ 
rent layer in the case with cr = 100 and domain size 
Lx X Lz = SOOdi X 194:di {Ly = SOOdi for the 3D sim¬ 
ulation) from runs 2D-7 and 3D-7. Panel (a) shows the 
color-coded current density from the 2D simulation and 
Panel (b) shows a 2 D cut of the current density and 
a 3D isosurface of plasma density colored by the cur¬ 
rent density from the 3D simulation at ujpet = 175 and 
ujpet = 375, respectively. Starting from the initial per¬ 
turbation, the current sheet gradually narrows as the 
current density is concentrated in the central region. 
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In the 2D simulation, the extended thin current sheet 
breaks into a number of fast-moving secondary plasmoids 
{ujpet ^ 225) due to the secondary tearing instability. 
The plasmoids coalesce and eventually merge into a sin¬ 
gle island at the edge of the simulation domain similar t o 
the nonrelativistic case (Daughton fc Karimabadi]|2QQ7). 

[d 


In the 3D simulation, as the intrinsic guide lieL 
ciated with the force-free current layer is expelled from 


the central region, the kink instability (Daughton 


1999) 


develops and interacts with the teying m ode, leading to 
a turbulent evolution (Yin et al. 2008). However, de¬ 
spite the strong 3D effects that modify the current layer, 
small-scale flux-rope-like structures with intense current 
density develop repeatedly as a result of the secondary 
tearing instability. 

Although the plasma dynamics in the 2 D and 3D simu¬ 
lations appears quite different, the energy conversion and 
particle energization are very similar. Figure 2 (a) shows 
the evolution of magnetic energy Eb^ electric field energy 
Ee^ kinetic energy and energy carried by relativistic 
particles with 7 > 4 from the 2D and 3D simulations (2D- 
7 and 3D-7). Note in both of these simulations, the total 
energy is conserved to within 10 “^ of the initial value. 
The evolutions of different forms of energies between 2D 
and 3D simulations are very similar. In both the 2D 
and 3D simulations, about 25% of the magnetic energy 
is converted into plasma kinetic energy, most of which 
is carried by relativistic particles. Figure 2 (b) shows 
the time-integrated energy conversion from magnetic en¬ 


ergy into plasma energy in the simulation dt f dVJ • E 
and its contribution from parallel and perpendicular elec¬ 
tric field terms Jy • Ey and • E^, respectively. Here 
f dV = f dxdydz. The difference in energy conversion 
between the 2D and 3D simulations can be as large as a 
factor of two at ujpet = 300, but at the end of the simula¬ 
tions both cases have converted about the same amount 
of magnetic energy. This shows that the kink instabil¬ 
ity that may modify the magnetic field does not signif¬ 
icantly change the overall energy conversion. While the 
energy conversion through parallel electric field is impor¬ 
tant when the thin current layer initially develops, most 
of the energy conversion is due to perpendicular electric 
fields induced by relativistic flows as the system is dom¬ 
inated by secondary plasmoids/flux ropes. This analysis 
has been done in all the cases and summarized in Table 
1 , which shows that in most of cases, the perpendicular 
electric field plays a dominant role in converting mag¬ 
netic energy into plasma kinetic energy. This can also 
be seen in Figure 3, which shows the color-coded inten¬ 
sities of J • E, Jx • Ex and Jy • Ey from the 2 D and 3D 
simulations at ujpet = 175 and ujpet = 375, respectively. 
Figure 2 (c) compares the energy spectra from the 2 D 
and 3D simulations at various times. The most striking 
feature is that a hard power-law spectrum / oc (7 — 1 )“^ 
with a spectral index p ^ 1.35 forms in both 2D and 
3D runs. Although a fraction of particles are acceler¬ 
ated in the early phase when the parallel electric field 
is important, most of the particles in the power-law dis¬ 
tribution are accelerated when the system is dominated 
by plasmoids/flux ropes. As we will discuss below, the 
formation of power law is closely related to the motional 
electric field induced by the fast moving plasmoids. In 
the subpanel, the energy spectrum for all particles in the 


3D simulation at ujpet = 700 is shown by the red line. 
The low-energy portion can be fitted by a Maxwellian 
distribution (black) and the nonthermal part resembles 
a power-law distribution (blue) starting at 7 ^ 2 with 
an exponential cut-off for 7 > 100. The nonthermal part 
contains ^ 25% of particles and ^ 95% of the kinetic 
energy. The maximum particle energy of the system can 
be predicted approximately using the reconnecting elec¬ 
tric field meC^{'ymax “ 1 ) = / |^^rec|cdt until the gyro- 
radius is comparable to the system size (see also Figure 
6 b). Although we observe a strong kink instability in 
the 3D simulations, the energy conversion and particle 
energy spectra are remarkably similar to the 2D results, 
indicating the 3D effects are not crucial for the particle 
acceleration. The fast acceleration is distinct from that 
of nonrelativistic magnetic reconnection, where particles 


are at most accelerated to mildly relativistic energy (o 


__ _ _ __ (e ,, 

Fu et al.1|2QQ6| [Drake et al.||2QQ6{ [Pritchett]|2Q06| ]Oka[ 
et al.[[2010[). The nonthermal-dominated distribution in 


the simulations is also quite different from distribution s 
in the relativistic shock regions (e.g., Spitkovsky 2008), 
where the particles are heated at the shock front arid form 
an extended thermal distribution containing most of the 
dissipated energy. The power-law spectral index p ^ 1 
from relativistic reconnection is significantly harder than 
the limit p ^ 2 predicted by nonrela tivistic and relativis¬ 


tic shock acceleration theorie s (e.g., Blandford & Eichler 
1987 Achterberg et al.[[ 2001 [ ) 


3.2. Particle Acceleration 


We now discuss the details of particle acceleration. We 
will first present some analysis of particle trajectories 
to show the acceleration mechanism. Then the domi¬ 
nant acceleration mechanism is distinguished by track¬ 
ing all the particles and calculating the energy gain us¬ 
ing the guiding-center drift approximation. The results 
demonstrate that the dominant acceleration mechanism 
is a first-order Fermi acceleration through curvature drift 
motion along the motional electric field induced by the 
relativistic reconnection flows. We calculate the acceler¬ 
ation rate a = Ae/{eAt) and its time integral for cases 
with cr = 6 — 400, where Ae is the averaged energy gain 
for particles of energy e over a period At. Finally, we 
summarize the character of the energy spectra. These 
main results will be discussed and interpreted in detail 
in Section 4, where we present the acceleration model. 

Figure 4 and Figure 5 present the trajectory analysis 
for the motions of accelerated particles in the 2D case 
with a = 100 and x Lz = 600di x 388di. These 
particles are selected to show the common acceleration 
pattern of accelerated particles, which is consistent with 
the results of statistical analysis for the acceleration of 
particles in Figure 6 . The first three panels of Figure 4 
show (a) the trajectory of a representative particle close 
to the central sheet between ujpet = 30 - 300 together 
with E|| at ujpet = 180, (b) the trajectory of the same 
particle between (jp^t = 310 - 510 together with Ey at 
ujpet = 400, and (c) the trajectory of the particle between 
(jpet = 510 - 720 together with Ey at ujpet = 640, respec¬ 
tively. The starting and ending locations of the particle 
are labeled by ‘-h’ and ‘x’ signs, respectively. Note that 
the field is highly variable in time and the location of 
the particle at the same time step as the field contour 
is drawn by the sign. The two bottom panels show 
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Fig. 1.— Evolution of 2D and 3D simulations with a = 100 and domain size LxX Lz = SOOd^ x 194d^ {Ly = SOOd^ for 3D); (a) Color-coded 
current density from the 2D simulation at (jOpet = 175 and oopet = 375, respectively; (b) 2D cut of current density and a 3D isosurface of 
the plasma density colored by the current density at oupet = 175 and oupet = 375, respectively. 


the evolution of the particle energy as a function of time 
(d) and energy as a function of the x position (e), re¬ 
spectively. Each period corresponding to that in (a)-(c) 
is labeled by the same color. The green curve represents 
the energy gain in the parallel electric field integrated 
from t = 0. Initially the particle is close to the central 
layer and gains energy by the parallel electric field. It is 
then strongly accelerated by perpendicular electric field 
when the reconnection region breaks into multiple islands 
and the electric field is mostly the motional electric field 
E = -V X B/c generated by relativistic plasma outflows. 
The figure also shows that the acceleration by resem¬ 
bles a Fermi process by bouncing back and forth within 
a magnetic island. 

Figure 5 presents another view of the particle accel¬ 
eration physics. It is similar to Figure 4, but the field 
contours show the outflow speed to highlight the role of 
Vx in the particle’s energization. This clearly illustrates a 
relativistic first-order Fermi process by bouncing in out¬ 
flow regions of the reconnection layer. Note the energy 
gain from the parallel electric field for this sample par¬ 
ticle is negligible since it entered the reconnection layer 
longer after the development of multiple plasmoids. 

In Figure 6, we present more analysis for the mecha¬ 
nism of particle acceleration. Panel (a) shows the energy 
as a function of the x-position of four accelerated parti¬ 
cles. Similar to Figure 5, the electrons gain energy by 
bouncing back and forth within the reconnection layer. 


We have analyzed trajectories of a large number of par¬ 
ticles and found the energy gain for each cycle is Ae ^ 5, 
which demonstrates that the acceleration mechanism is a 
first-order F ermi process ( Drake et al.|2QQ6 2010 Kowal 


et al. 2011). Panel (b) shows the maximum particle en¬ 


ergy in the system as a function of time. This is plot 
ted using different count level from the 1-particle level 
to the 1000-particle level. Also plotted is the estimated 
maximum energy resulting from the reconnecting electric 
field by assuming particles moving along the electric field 
at the speed of light f \qErec\cdt. This shows that the 
maximum possible energy occurs for a small number of 
particles that continuously sample the reconnection elec¬ 
tric field rrieC^^fmax = f iQ-E^recl^dt. At late time, as the 
particle gyroradius becomes large and comparable to the 
system size, the maximum energy saturates. To show 
the Fermi process more rigorously, we have tracked the 
energy change for all the particles in the simulation and 
the relative contributions arising from the parallel elec 
trie field (rUeC^Ay = f qvjjEndt) and curvature drift ac 


celeration {rrieC^Aj = f qVcurv similar to (Dahlin 

et al.|[2Ql4 ), where Veurv = 7'^y (b x (b • V)h)/Qce, '^|| is 


the particle velocity parallel to the magnetic field, and 
b = 'B/\B\. Panel (c) shows the averaged energy gain 
and the contribution from parallel electric field and cur¬ 
vature drift acceleration over an interval of 2buj~^ as a 


function of energy starting at ujpet = 350. The energy 
gain follows Ae ^ as^ confirming the first-order Fermi 
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CDpet 



copet 



y-1 


Fig. 2.— Plasma energetics in 2D and 3D simulations with 
cr = 100 and domain size Lx X Lz = SOOdi x 194di {Ly = 300di for 
3D); (a) Evolution of magnetic energy Eb^ electric field energy Ee^ 
plasma kinetic energy Ej^ and energy carried by relativistic parti¬ 
cles with Lorentz factor 7 > 4; (b) Energy conversion from mag¬ 
netic energy into plasma energy integrated over time dt J dV J • E 
and its contribution from parallel and perpendicular electric field 
JII • E|| and (c) Evolution of particle energy spectra from 

21) and 3D simulations. Subpanel: Energy spectrum from the 3D 
simulations at Upet = 700. The low energy is fitted with a thermal 
distribution and rest of the distribution is a nonthermal power law 
with an exponential cutoff. 


process identified from particle trajectories. The energy 
gain from the parallel motion depends weakly on energy, 
whereas the energy gain from the curvature drift accel¬ 
eration is roughly proportional to energy. In the early 
phase, the parallel electric field is strong but only accel¬ 
erates a small portion of particles, and the curvature drift 
dominates the acceleration starting at about ujpet = 250. 
The contribution from the gradient drift was also evalu¬ 
ated and found to be negligible in comparison. Panel (d) 
shows a =< As > /{sAt) measured directly from the 
energy gain of the particles in the perpendicular elec¬ 
tric field (rUeC^Ay = j qw^ * and estimated from 

the expression for the curvature drift acceleration. The 
close agreement demonstrates that curvature drift term 
dominates the particle energization. 

For higher a and larger domains, the acceleration is 
stronger and reconnection is sustained over a longer du¬ 
ration. In Figure 7(a), we present the energy spectra at 
the end of simulation for a number of cases with different 
(7 and system size x = 600d^ x SSSdi. A summary 
for the spectral index can be found in Table 1. In Fig¬ 
ure 7(b), a summary for the measured spectral index for 
the power-law ranges of all the 2D runs shows that the 
spectrum is harder for higher a and larger domain sizes, 
and approaches the limit p = 1. Note that the spec¬ 
tral indexes a ppear systematically harder than in other 


recent papers (Sironi fc Spitkovsky 2014 Melzani et al 


2014b Werner et al.||2014p . However, the energy spectra 


in these studies are plotted using total relativistic en¬ 
ergy ymc^ and here we use kinetic energy (7 — l)mc^. 
Using total relativistic energy in the energy spectra sig¬ 
nificantly distorts the spectral index in the energy range 
of 0 < 7 — 1 < 10 , which may alter the interpretation 
of the results ([Sironi fc Spitkovsky 2014 Melzani et al. 
2014b Werner et al.||2014ip| 


3.3. Reconnection rate and relativistic flows 

Figure 8 (a) shows the time-dependent reconnection 
rates normalized using the initial asymptotic magnetic 
field Bq in 2D and 3D simulation with a = 100 (Run 
2D-7 and 3D-7). The 2 D reconnection rate is computed 
from 


_ Epee 1 d^Jj 

^ ^ 50^40 

where 7 /) = max(A^) — min(A^) along the central layer 
z = 0, Ay is the vector potential along the y direction 


<> represents a time average over Stujpe = 25 (Liu et al. 


2014), Vao = va/\/1 + {va/cY = is the 

relativistic Alfven speed in the cold-plasma limit. Here 
va = non-relativistic Alfven 

speed based on 5o. The 3D reconnection rate is esti¬ 
mated by using the mixing of plasma across the sepa- 
ratrix surfaces (Daughton et al. 2014). The rate in the 


^ In fact, our simulation results show that the “— 1 ” spectra can 
be obtained as long as the magnetic energy dominates over the 
initial plasma kinetic energy 87 rnfcTo/B^ = /3 <C 1. An example 
can be seen in the Appendix (Figure 13), which robustly shows the 
p = 1 spectrum can be obtained when cr = 25 and A:To = O.Olmc^. 
The same spectrum gives a “p ~ 2” slope when it is plotted as a 
function of 7 , wh ich may explain the different conclusions reporteci 
by other papers JS ironi fc Spitkovskvll2014| |Melzani et al.||2014b| 
[Werner et ai.|201^ 
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Fig. 3.— Color-coded intensity of energy conversion rate J • E normalized using nomeC^(^pe and contributions from J_l • E_l and Jy • Ey 
for the 2D and 3D simulations with a = 100 at oopet = 175 and oopet = 375, respectively. In the early stage the conversion by parallel 
electric field is important and the perpendicular electric field plays a dominant role when multiple-plasmoids (fiux ropes in 3D) develop 
due to the secondary tearing instability. 


2D simulation is quite variable but the range is within a 
factor of two times of the 3D results, meaning that the 
2D and 3D simulations give roughly the same reconnec¬ 
tion rate. Figure 8 (b) shows the peak reconnection rate 
for a number of 2D cases with a from 0.25 to 1600 and 
box size 1200di x 776d^. The rate is observed to increase 
with a from Epee ^ 0.035o for cr = 1 to Epee ^ 0.245o 
for a = 1600. It shows that the peak reconnection field 
increases with a and starts to saturate around a = 1000. 
For low-cr cases with cr < 1, the reconnecting electric 
field is consistent with previous work for nonrela tivistic 
reconnection (e.g., Daughton & Karimabadi|2007). More 
detailed analyses have shown that for high-cr cases, the 
reconnection rate normalized using the magnetic field Bu 
upstrea m of the diffusion region Ep^dBu is close to 1 for 
cr > 100 dLiu et al.|[M^ . 

In Figure 9 (a) and (b) we plot the maximum flow 
velocity in the x direction (outflow direction) and the 
corresponding Lorentz factor F^,. The 2D results are 
represented by blue symbols and the 3D results are in 
red symbols, respectively. Although we have only used a 


small simulation domain that may be affected by counter¬ 
streaming particles, a relativistic outflow still develops 
with Vx of a few. In Figure 9 (c) and (d) we plot the max¬ 
imum flow velocity in the z direction (inflow direction) 
and the corresponding Lorentz factor respectively. 
Interestingly, the inflow speed can also be relativistic for 
high-cr cases. Detail ed analysis for t he diffusion region 
has been discussed in Liu et al. (2015), which shows that 
the inflow speed can be predicted by a model based on 
the Lorentz contraction of the plasma passing through 
the diffusion region. 

The enhanced reconnection rate and development of 
relativistic inflow/outflow s tructures are in contrast t o 
the results reported earlier (Sironi & Spitkovsky 2014), 
where the outflow can only be mildly relativisti c and the 
inflow speed remains nonrelativistic. Note that |Liu et al. 
(2015) has also reported the development of relativis 


tic inflow for both Harris and force-free current sheets, 
indicating that this property of relativistic magnetic re¬ 
connection does not strongly depend on the initial setup. 
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Fig. 4. — Panels (a)-(c) show a particle trajectory in the x-z plane together with the color-coded electric field (a) , (b) Ey^ and (c) 

Ey. Panels (d) and (e) show the particle energy as a function of time and energy as a function of the x position, respectively. In (d) and 
(e), curves with different colors represent the energy evolution during time periods in (a)-(c). The green curve shows the integrated energy 
gain from the parallel electric field. 


3.4. 3D Dynamics 

In our three-dimensional simulation, we also find 
strong bulk ^ 4 can develop in the system, meaning 
the development of relativistic flows is not strongly influ¬ 
enced by 3D effects. Figure 10 shows the power spectrum 
of magnetic fluctuations with wave numbers perpendic¬ 
ular to the y direction and a volume rendering of the 
current density in the 3D simulation with a = 100 at 
ujpei^ = 708. The power spectrum shows a clear iner¬ 
tial range with a slope of “—2” and steeper slope for 
higher wave numbers k^di ^ 1. As we have discussed, 
the 3D simulation allows the development and interac¬ 
tion of secondary tearing instability and kink instability, 
leading to a turbulent magnetic field in the reconnection 
layer. For the whole simulation, the range of scales for 


the 2D magnetic islands is similar to the observed 3D 
flux ropes. The maximum energy in both 2D and 3D 
agrees well with the time integral of energy gain from 
reconnecting electric fie ld. This is in contrast to the ear¬ 
lier k inetic simulations ( Zenitani fc Hoshino|2QQ5~ " 


2007 


2008). The energy distributions reported in this paper 


are remarkably similar in 2D and 3D, suggesting that the 
underlying Fermi acceleration is rather robust and does 
not depend on the existence of well-defined magnetic is¬ 
lands. 


4. A SIMPLE MODEL 

It is often argued that some loss mechanism is needed 
to fo r m a power-law distribution (|Zenitani fc Hoshino 


2QQT| [Drake et al^lMol |HQshino||2012[ ). How- 
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Fig. 5. — Panels (a)-(d) show a particle trajectory in the x-z plane together with the fluid velocity in the x direction Vx- Panels (e) and 
(f) show the particle energy as a function of time and energy as a function of the x position, respectively. Different colored curves represent 
the energy evolution during time periods in (a)-(d), showing that the particle gains energy by bouncing in the relativistic ffow generated 
by reconnection. 


ever, the simulation results reported in this paper clearly 
show power-law distributions in a closed periodic sys¬ 
tem. Here we present a simple model to explain the 
power-law energy spectrum observed in our PIC simu¬ 
lations. The model is illustrated by Figure 11(a). As 
reconnection proceeds, cold plasma in the upstream re¬ 
gion advects into the acceleration zone at a constant ve¬ 
locity that is determined by reconnection electric field 
Vin = cE^ec X B/H^. The process lasts r ^ 
where is the size of the simulation box along the 2 : 
direction. In the acceleration region, our analysis has 
shown that a first-order Fermi process dominates the 
energy gain during reconnection. We solve the energy- 
continuity equation for the energy distribution function 


/(5,t) within the acceleration region 


df d 



( 1 ) 


with acceleration dejdt = ae^ where 5 = meC^{j — l)/kT 
is the normalized kinetic energy and a is the constant 
acceleration rate in the first-order Fermi process. We 
assume that the initial distribution within the layer /o 
is Maxwellian with initial temperature kT < rUeC^, such 
that 


/o(X 7(72 - 1)^/2 exp(-e) (2) 
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^ 1400 





Fig. 6. — (a) Energy as a function of ar-position of four accelerated particles; (b) The maximum energy of particles in the system as a 
function of time from the 1-particle count level to the 1000-particle count level. The red dashed line shows the maximum energy estimated 
for a particle moving in the direction of the reconnecting electric field at the speed of light nieC^^max = f \qErec\cdt] (c) Averaged energy 
gain and contributions from parallel electric fields and curvature drift acceleration over a time interval of 2bujpe as a function of particle 
energy starting at (jUpet = 350;(d) a =< As > /{sAt) from energy gain in perpendicular electric field and by curvature drift acceleration, 
and from the Equation (6) using the averaged fiow speed and island size. 


Run 

a 

system size 

A 

P 

'ymax 

Ekin^^ 

(J-E)±% 


2D-1 

6 

SOOdi X WMi 

Odi 

2.2 

45 

23% 

83% 

0.4 

2D-2 

6 

600di X SSSdi 

Odi 

2.0 

56 

32% 

92% 

0.5 

2D-3 

6 

1200di X 776di 

6di 

1.7 

79 

34% 

93% 

0.7 

2D-4 

25 

SOOdi X lOAdi 

Odi 

1.6 

195 

28% 

85% 

1.1 

2D-5 

25 

OOOdi X SSSdi 

Odi 

1.3 

339 

37% 

90% 

1.6 

2D-6 

25 

1200di X 776di 

Odi 

1.2 

617 

42% 

90% 

2.0 

2D-7 

100 

SOOdi X lOAdi 

Odi 

1.35 

650 

29% 

73% 

2.0 

3D-7 

100 

SOOdi X lOMi X SOOdi 

Odi 

1.35 

617 

28% 

71% 

N/A 

2D-8 

100 

OOOdi X SSSdi 

Odi 

1.25 

1148 

40% 

78% 

3.1 

2D-9 

100 

1200di X 776di 

6di 

1.15 

1862 

45% 

94% 

4.3 

2D-10 

400 

SOOdi X 194di 

12di 

1.25 

1514 

20 % 

54% 

3.0 

2D-11 

400 

OOOdi X SSSdi 

12di 

1.15 

3715 

31% 

75% 

4.8 

2D-12 

400 

1200di X 776di 

12di 

1.1 

5495 

36% 

86 % 

6.5 

2D-13 

1600 

SOOdi X 194:di 

24di 

1.2 

2812 

13% 

45% 

N/A 

2D-14 

1600 

OOOdi X SSSdi 

24:di 

1.1 

7913 

21 % 

53% 

N/A 

2D-15 

1600 

1200di X 776di 

24di 

1.05 

11220 

30% 

66 % 

N/A 


TABLE 1 

List of simulation runs with ct ^ 6 . The spectral index p, the maximum energy (100- particle level) at the end of the 

SIMULATION 7macc, THE PERCENTAGE OF MAGNETIC ENERGY THAT IS CONVERTED INTO KINETIC ENERGY THE CONVERSION OF 

MAGNETIC ENERGY CAUSED BY PERPENDICULAR ELECTRIC FIELD ( J • E) ±_ AND CXTinj ESTIMATED BY TRACKING PARTICLES IN THE SYSTEM. 


For simplicity, we consider the lowest order (nonrela- 
tivistic) term in this expansion and normalize /o = 
^^V^exp(— 5 ) by the number of particles A^o within the 
initial layer. The distribution after time t is: 


/(e,t) = ^Vee ^“*/7xp(-ee “*), (3) 


which remains a thermal distribution with a tempera¬ 


ture consistent with that obtained by Drake et al. 


(|2010|). However, since upstream particles enter contin¬ 
uously into the acceleration region, the number of par¬ 
ticles in the acceleration zone increases with time. We 
consider a particle distribution fi^j = v^exp(—g) 

with number of particles Ni^j oc VinTinj injected from 
upstream, where Tinj is the time scale for particle injec- 
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Fig. 7.— Energy spectra at the end of simulations for a series of 
2D runs with system size Lx x Lz = 600di x 388di and different 
(7 from 6 to 1600. (b) Spectra index for all 2D simulations with a 
from 6 to 1600. (c) Time integrated arinj for cases with a = 6-400 
and different system sizes. 


tion. To highlight the key role that time-dependent injec¬ 
tion plays in setting up the power-law, we first consider a 
quick heuristic derivation of the main result. To proceed, 
we split finj into N groups and release the group into 
the acceleration region at time t = jAt. Since each group 
will satisfy the Equation^ for a different initial time, af¬ 
ter we have injected the final group at t = Tinj^ the total 
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Fig. 8 .— (a) Time-dependent 2D and 3D reconnection electric 
field normalized by the initial magnetic field Erec/Bo. (b) Nor¬ 
malized peak electric field EreclBo as a function of a in 2D simu¬ 
lations. 


distribution (in the limit N ^ oo) is 

/(£, t) - {-ee-^*)dt (4) 

y'^Enj Jo 

Ninj jerf(52/^) — erf(52/^e“^'^"’^-^/2) 

^ A 

In the limit of ar 1, this gives the relation f oc 1/e 
in the energy range 1 < s < Figure [T^b) shows 

(4) for different arinj • A power-law spectrum with p = 1 
emerges as ar^nj increases ar^nj > 1. Note that for a 
closed system, since the averaged magnetic energy per 
particle is only ameC^/d and the energy in each energy 
bin is constant, the maximum energy the power law can 
only extend to ^rnax ^ cr/4. 

Next, in order to treat the problem more rigorously, 
and include the inffuence of particle escape, we consider 
the more complete equation 

E + E (^f \ = El _ _L_ 

dt de\dt^J Tinj Tesc' 


( 5 ) 
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Fig. 9.— (a) The maximum flow velocity in the x direction Vx as a function of cr; (b) The maximum flow Lorentz factor in the x direction 
Fx = 1/(1 — as a function of cr; (c) The maximum flow velocity in the 2 ; direction I 4 as a function of cr; (d) The maximum flow 

Lorentz factor in the 2 ; direction Fz = 1/(1 — jc^) as a function of cr. 


where Tesc is the escape time for particles. For the initial 
current-layer distribution /o and injected particle distri¬ 
bution finj considered above, the solution can be written 
as 


f{e,t) = —^y/ee 

V 71 


-(3/2+/3)at ^ 


K.F{-ee 


(6) 


+ 


2N- ■ 


^{aTinj)e'^+^ 


[r(3/2+/3)(£e “*) - r(3/2+;3)(e)] , 


where P = l/{aTesc) and Ts{x) is the incomplete Gamma 
function. The first term accounts for particles initially in 
the acceleration region while the second term describes 
the evolution of injected particles. In the limit of no in¬ 
jection or escape {resc 00 and Tinj oc), the first 
term in remains a thermal distribution the same as 
(3). However, as reconnection proceeds new particles en¬ 
ter continuously into the acceleration region and due to 
the periodic boundary conditions there is no particle es¬ 
cape. Thus considering the case Tesc ^ and assuming 
Nq Ninj , at the time t = Tinj when reconnection satu¬ 
rates the second term in (6) simplifies to (4). Thus in the 
limit Nq Ninj the first term in (6) should be retained 
and the power-law produced is sub-thermal relative to 
this population. While it is straightforward to obtain the 
relativistic corrections arising from the injected distribu¬ 
tion (2), we emphasize that these terms do not alter the 
spectral index. This solution explains results from our 
simulations, and also appears to explain the results from 
several recent papers, which obtained power-law distri¬ 


butions by subtracti ng the initial hot plasma component 


in the current layer (ISironi & Spitkovsky 2014 

Melzani 

et al. 

2014b 

Werner et al.| 

|20I4|). In particular. 

IVIelzani 


cle distribution initially in the current layer and reported 
it as a heated Maxwellian distribution. 

In order to estimate the acceleration rate a, the en¬ 
ergy change of each particle can be approxima ted by a 
relativistic collision formula (e.g., Longair|[l994 ) 


Ae = r^(l + 


2Vv^ 






(7) 


where V is the outflow speed, Ty = 1/(1 — and 

Vx is the particle velocity in the x direction. The time be¬ 
tween two collisions is about Lisjvx^ where Lis is the typ¬ 
ical size of the magnetic islands (or flux ropes in 3D). As¬ 
suming that relativistic particles have a nearly isotropic 
distribution Vx ^ c/2, then 


Ae c(rf.(l + ^ + ^)-l) 
sAt ‘^Lis 


(8) 


Using this expression, we measure the averaged V 
and Lis from the simulations and estimate the time- 
dependent acceleration rate a(t). An example is shown 
in Figure 6 (d). This agrees reasonably well with that 
obtained from perpendicular acceleration and curvature 
drift acceleration. Figure 7(c) shows the time-integrated 
value of aTinj = a{t)dt for various simulations with 
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Fig. 10. — The evidence for turbulence in the 3D simulation. Left: power spectrum of magnetic fluctuations with wave numbers 
perpendicular to the y direction. Right: volume rendering of the current density J/Jq in the 3D simulation at uopet = 708. 


cr = 6 — 400. For cases with aTinj > 1, a hard power- 
law distribution with spectral index p ^ 1 forms. For 
higher cr and larger system size, the magnitude of ar^nj 
increases approximately as oc cr^/^. 

Better agreement between the simple model and the 
PIC simulations can be reached by considering the time- 
dependent acceleration rate a{t). As the magnetic recon¬ 
nection saturates, the acceleration rate decreases. Figure 
11 (c) shows the solution that uses the time-dependent 
acceleration rate a{t) in Figure 6 ( d) using a stochas tic in¬ 
tegration technique described by|Guo et al.| (2010). The 
final spectral index is about p = 1.25, similar to that 
from the PIC simulation shown in Figure 7(a). 

5. IMPLICATIONS 

We discuss the implication of the above conclusions for 
understanding the role of magnetic reconnection in mag¬ 
netically dominated astrophysical systems. Based on the 
current understanding of magnetic reconnection, mul¬ 
tiple X-line reconnection develops when the secondary 
tearing instability is active in large-scale collisionless 
plasma system. This process may also be important 
when a hierarchy of collisional plasmoids (Loureiro et al. 


20071 Bhattacharjee et al. 2009] jUzdensky et al.||2010) 

1 « 


develops kinetic scale cur rent layers that may trigger 


collisionless reco nnection (Daughton et al. 2009 Ji & 


Daughton 2011|). Therefore the collisionless reconnec¬ 
tion process discussed here is relevant to a range of high- 
energ y astrophysical problems below (seejJi & Daughton 

■y of astrophysical 


2011 [ for a comprehensive summary 
problems with relevant physics). 

5.1. Pulsar Wind Nebulae 

In PWN models, magnetic reconnection has been 


ergy in Poynting-flux dominated flows (|Coroniti| 

llpgo 

Lyubarsky & Kirk||20011 Kirk & Skiaeraasen||2UU3t 

Forth 

et al. 

2013 

) and accelerating particles to high energies 
). In PWNe, the emission flux usually has 

(|Kirk 

2DD3 


requires an electron energy distribution dN/dj oc 7 “^ 
with p = 1 — 1.6 {p = 2ay + 1 ), too hard to be ex - 
plained by diffusive shock acceleration (|Atoyan||1999). 
The recently detected > 100-MeV Crab flares have pho- 
ton energies well above the usually employed upper limit 
for synchrotron emissions, challenging the traditional ac- 
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Fig. 11.— (a) Illustration of the acceleration model for the forma¬ 
tion of power-law distributions; (b) Analytical results for different 
Ci'Tinj obtained from Eq. (4). (c) The solution of Eq. (1) using 
time-dependent a{t) from Figure 6 (d). 


celeration theory ([Abdo e t al. 2011 Tavani et al. 2011 
Biihler & Blandford|2014|). T'liere are two mam possibil- 
me photi 


ities tor explaining tne photon energies, ( 1 ) a relativistic 


Doppler boosting of the emitting region (Clausen-Brown 


& Lyutikov 2012 ) and/or ( 2 ) a strong particle acceler- 


ation m a nonideal electric field where E > where 
is the magnetic field perpendicular to pa rticle veloc¬ 
ity (iCerutti et al.j|2013t |Lyutikov et al.||2oi4 ). 

Tnese observations suggest that relativistic magnetic 
reconnection may occur in the Crab nebula. The power 
law index revealed in this study is p = 1 — 2 , consis 


tent w ith the inferred spectra in the radio range (Atoyan 
1999[) and in high-e nergy during the Crab 7 —ray flares 


Vl'avani et al. 2011). Explaining these observations re- 
quires a fast and efficient dissipation mechanism that 
converts a substanti al fraction of niagnet ic energy into 
relativistic particles ( Lyutikov et al.|[2014 ). In the Crab 


pulsar, magnetic reconnection is estimated to be in the 
plasmoid dominated regime and can dissipa te a nontriv¬ 


ial fraction of th e pulsar spin-down power (Uzdensky & 


Spitkovsky|2014). Our simulations have shown that for a 


magnetically dominated reconnection layer with a ^ 1 
magnetic reconnection rate is greatly enhanced by about 
one order of m agnitude compar ed to the nonrelativistic 
limit (see also Liu et al. 2015) and a large fraction of 
magnetic energy in the system is converted into nonther- 
mal energy distribution, suggesting an efficient magnetic 
dissipation and strong nonthermal radiation processes in 
the Crab wind. The maximum particle energy increases 
linearly and can be well predicted by assuming particles 
moving along the reconnecting electric field at the speed 
of light. There are also relativistic inflow and outflow 
structures {Tmax ^ 10 ) associated with reconnection, 
which may boost the emission photo n energy and help 


to explain the observed Crab flares (Clausen-Brown & 
Lyutiko^| 2 Q 12 ). It is interesting to note that the recon- 


nection acceleration may also explain the pulsed 7 -ray 
emission, although observations at high er energies is re¬ 


quire d to further constrain the model (Mochol & Petri 


2015). 


5.2. AGN Jets 

In AGN jets, a number of 7 -ray sources have flat radio 
spectra with indices around = 0 , meaning the el ectron 
energy distribution index may be close to p = 1 (Abdo 
et al.||2010 Hayashida et al. 2015). Several blazars have 
shown extremely fast va riability in TeV range on the or¬ 
der of seve ral minutes (Aharonian et al. 2007 Albert 
et al. 2007). Hard power laws p ^ 1 in T'eV range have 
been inferred after removing effect of the extragalactic 
background light using variou s models (Aharonian et al. 
20061 Krennrich et al. 2008). For GeV-TeV flat spec- 
trum rad io quasars (F)SKCj ) , high radiation efficiency is 
reported (Zhang et al.|2013) and the electron (Je which is 
measured as magnetic energy power to the elect ron en¬ 
ergy power is very high up to the order of 100 (Zhang 
et al.JboMD. 


Explaining the fast variability requires the relativis¬ 
tic beaming effect p ossibly arising from relativistic re¬ 
conne ction outflows (Giannios et al. 2009 Deng et al. 
2015). Our kinetic simulations have shown that the 
Lorentz factor of the maximum outflow speed T^, ^ 10 
for a 1000. The simulation results and theoretical 
model predict hard particle energy distribution consis- 
tent with the hard radio spect ra observed in some AGNs 
(Romanova & Lovelace 1992). Recent advanced AGN 
emission models have interred that at least for some 
types of blazers, particularly FSRQ, strong particle ac¬ 
celeration and/or strong magnetic field is necessary to 
explain fast flares and a inferred fro m the model fittin g 
can significantly exceed unity a ^ 1 (Chen et al. 12014) 


Magnetic reconnection may offer an explanation for the 
simultaneous decrease of magnetic field and emission in¬ 
crease during the flare phase of blazar flares a nd is a 
promising s cenario for modeling AGN emissions (Zhang 


5.3. Gamma-ray bursts 

In gamma-ray bursts (GRBs), the traditional internal 
shock model of prompt emission is difficult to reconcile 
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with observations (see Zhang & Yan|2011 and references 
therein). Magnetic reconnection and associated particle 
acceleration have been proposed as a key process in GRB 
models such as ICMART mod el dZhang fc Yan||2Qll|) 
and r econnection-switch model (IVlcKinney Uzdensky 
2012). The efficient magnetic dissipation and particle 
acceleration during reconnection may be impor tant to 


understand the emission mechanism in GRBs ([Kumar 
1999| Spruit et al. 2001 Drenkhahn & Spruit 2002^7 Gru¬ 


ber et al.| (|2014p have shown a series of features in GRB 
prompt emission that are not consistent with the simple 
synchrotron shock model. For example, the hard low- 
energy spectra, where the particle energy spectral index 
is close to p = 1 assuming synchr otron radiation (Ghis- 
ellini et al. 2000 Preece et al.| 2002 ) and the thermal emis¬ 
sion component predicted in the ffieba ll-internal-shock 
model has been rarely seen in GRBs (Zhang & Pe’er 
200 ^|Abdo et al.|2009|). 

From our simulation results and analytical model, the 
particle energy spectral index is close to p = 1 , con¬ 
sistent with low-energy photon spectra observed in most 
GRB s ( Band et al.|1993 Preece et al .|2000 Gruber et al. 
2014). The acceleration in reconnection layers is much 
taster than the radiation cooling and can maintain the 
hard spectrum. Using PIC simulation, Spitkovsky (2008) 
found that in the downstream region of highly relativis- 
tic shocks the number of particles in the nonthermal tail 
is ~ 1 % of the entire downstream population, and they 
carry 10 % of the kinetic energy in the downstream re¬ 
gion. In our simulations of relativistic reconnection, the 
number of nonthermal relativistic particles is ^ 25% of 
the total number particles in the simulation and they 
carry ^ 95% of kinetic energy in the system, meaning 
relativistic reconnection is much more efficient in pro¬ 
ducing nonthermal relativistic particles. This efficient 
conversion from magnetic energy into kinetic energy of 
nonthermal particle s may help solve the efficiency prob - 
lem in GRB models ( Zhang et al.[2QQ7 ; Deng et al.|2Q15 ). 


5.4. Nonrelativistic reconnection sites 

While the primary focus of this paper is relativistic 
magnetic reconnection, the physics of Fermi acceleration 
and the formation of power-law distribution is also appli- 
cable to the nonrelativist i c regi mes previously discussed 


(Drake et al. 2006 2010 2013). Based on our analyti- 


cal model, the power-law distribution forms only when 
> 1 - The results in this paper demonstrate that 
this condition is more easily achieved in regimes with 
cr ^ 1 , but it may also occur with a < 1 in sufficiently 
large reconnection layers. In several preliminary simula¬ 
tions, we have observed the formation of similar power 
laws in nonrelativistic proton-electron plasma and will 
report elsewhere. 

X-ray observations of solar flares have shown strong 
particle acceleration and energy conversion during mag¬ 
netic reconnection and the particle distribution often 
takes power-law distributions, requiring a particle ac¬ 
celeration mecha nism that is dominated by nonther¬ 
mal acceleration (Krucker et al. 2008 [2010 Krucker & 

hi 


Battaglia 2014). As we have shown here, in magnet¬ 
ically dominated regimes, a large fraction of magnetic 
energy can be converted into particles in a power-law 
distribution. Similar process is likely to occur in so¬ 
lar flares, where the plasma [3 = ^T:nkTlB‘^ ^ 0.001 - 


0.01 (cr < 1 ). However, physics such as the influence 
of mijm^^ strong trapping at X-line region, and particle 
escape from the system need to be investigated further 
(Egedal & Daughton 2015, in preparation). 

6 . DISCUSSION AND CONCLUSION 

The dissipation of magnetic field and particle energiza¬ 
tion in the magnetically dominated systems is of strong 
interest in high energy astrophysics. In this study, we use 
2D and 3D fully kinetic simulations that resolve the full 
range of plasma physics to investigate the particle accel¬ 
eration and plasma dynamics during collisionless mag¬ 
netic reconnection in a pair plasma with magnetization 
parameter a varying from 0.25 to 1600. A force-free cur¬ 
rent layer, which does not require a hot plasma popu¬ 
lation in the current layer, is implemented as the initial 
condition. 

We find that the evolution of the current sheet and ac¬ 
celeration of particles has two stages. In the early stage, 
an extended reconnection region forms and generates a 
parallel electric field that accelerates particles in the cur¬ 
rent layer. As time proceeds, the layer breaks into mul¬ 
tiple plasmoids (flux ropes in 3D) due to the secondary 
tearing instability. The motional electric field in the re¬ 
connection layer strongly accelerates energetic particles 
via a first-order relativistic Fermi process leading to the 
conversion of most of the free energy in the system. A 
large fraction of the magnetic energy is quickly converted 
into the kinetic energy of nonthermal relativistic parti¬ 
cles (within a few light-crossing times) and the eventual 
energy spectra show a power law / oc (7 — 1 )“^, with 
the spectral index p decreasing with a and system size 
and approaching p = 1. The formation of the power- 
law distribution can be described by a simple model that 
includes both inflow and the Fermi acceleration. This 


(ISironi & Spitkovsky||20I4 

Melzani et al.| 

|20I4b 

|Werner 

et al. 2UI4p, which reported hard power-law d 

Listribu- 


tion inside the current layer. For the more realistic limit 
with both particle loss and injection, the spectral index 
p = 11 /{(yTesc): recovering the classical Fermi solution. 
If the escape is caused by convection out of the recon¬ 
nection region Tesc = Lx/Vx^ the spectral index should 
approach p = 1 when aTesc ^ 1 in the high-a regime. In 
preliminary 2D simulations using open boundary condi¬ 
tions, we have confirmed this trend and will report else¬ 
where. For the nonrelativistic limit, the reconnection 
needs to be sustained over a longer time to form a power 
law. 

We have also shown that in sufficiently high-a regimes 
the magnetic reconnection rate is enhanced and relativis¬ 
tic inflow and outflow structures develop. The scaling 
follows the prediction based on the Lorentz contraction 
of plasma passing through the diffusion region. Although 
3D magnetic turbulence is generated as a consequence of 
the growth of the secondary tearing instability and kink 
instability, the particle acceleration, energy release and 
reconnection rate in the 3D simulation are comparable 
to the corresponding 2D simulation. 

Our study has demonstrated that relativistic mag¬ 
netic reconnection is a highly efficient energy-dissipation 
mechanism in the magnetically dominated regimes. The 
plasma distribution in the reconnection layer features 
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power-law energy spectra, which may be important in 
understanding the nonthermal emissions from objects 
like pulsars, jets from black holes, and gamma-ray bursts. 
Both the inflow and outflow speeds approach the speed 
of light and have Lorentz factors of a few, which may ex¬ 
plain the fast variability and high luminosity observed in 
those high-energy astrophysical systems. These flndings 
on particle acceleration and plasma dynamics during rel¬ 
ativistic reconnection substantiate the important role of 
magnetic reconnection in high-energy astrophysical sys¬ 
tems. 
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APPENDIX: NUMERICAL CONVERGENCE 


The accuracy of particle-in-cell (PIC) kinetic simula¬ 
tions depends on a series of numerical parameters such 
as cell size, ti me step, and the number o f macro-particles 


per cell (e.g., 

Birdsall & Langdon 

1991). The numerical 

convergence of simulation results 

aas been rarely explic- 


itly checked when modeling astrophysical problems using 
PIC simulations, and often a small number of macro¬ 
particles are used. Here we examine the numerical con¬ 
vergence of our results on these numerical parameters 
using VPIC code for different initial temperatures from 
kTo = 0.01 to 0.36 rrieC^ . Our test case has a = 25 with 
box size x = 600d^ x SSSdi and simulation time 
ojpet = 3000. We And that numerical heating can become 
unacceptably large when a small number of particles per 
cell is used. In Table 2 we list the key parameters for the 
test. Although for most cases, the violation in energy 
conservation is small (Egrr/^tot within 1%), the numeri¬ 
cal heating can signiflcantly modify the particle distribu¬ 
tion since the initial kinetic energy is a small fraction of 
the total energy. Therefore to obtain trustworthy results 
that are numerically converged, the violation of energy 
conservation should be much less than the initial kinetic 
energy E’grr/^'/cO ^ 1- Eigure 13 shows several cases 
with kTo = 0.36meC^ with grid number 2048 x 2048, 
Courant number = 0.7, and different numbers of par¬ 
ticles per cell from 8 to 512. Eigure 14 shows several cases 
for kTo = O.OlmgC^ with grid number 4096 x 4096 and 
Cr = 0.9 but different numbers of particles per cell from 2 
to 512. Both flgures show that as the total energy change 
in the numerical simulations becomes smaller than the 
initial kinetic energy ^ 1, the numerical heat¬ 

ing has a negligible effect on the distribution function. 
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Run 

kTolnieC^ 

Grid numbers 

Time step (Gr) 

NPG 

Eerr ! EfQtal 

Eerr ! 

A-1 

0.36 

4096 X 4096 

0.9 

2 

10 % 

159% 

A-2 

0.36 

4096 X 4096 

0.9 

8 

2.4% 

38% 

A-3 

0.36 

4096 X 4096 

0.9 

32 

0.56% 

9% 

A-4 

0.36 

4096 X 4096 

0.9 

128 

0.084% 

1.3% 

A-5 

0.36 

2048 X 2048 

0.9 

8 

5% 

80% 

A-6 

0.36 

2048 X 2048 

0.9 

32 

1 .2% 

20 % 

A-7 

0.36 

2048 X 2048 

0.9 

128 

0.3% 

5% 

A-8 

0.36 

2048 X 2048 

0.7 

8 

1.9% 

30% 

A-9 

0.36 

2048 X 2048 

0.7 

32 

0.45% 

7% 

A-10 

0.36 

2048 X 2048 

0.7 

128 

0 .12% 

1.9% 

A-11 

0.36 

2048 X 2048 

0.7 

512 

0.04% 

0 .6% 

A-12 

0.36 

2048 X 2048 

0.5 

8 

0.75% 

12 % 

A-13 

0.36 

2048 X 2048 

0.5 

32 

0.19% 

3% 

A-14 

0.36 

2048 X 2048 

0.5 

128 

0.05% 

0 .8% 

B-1 

0.09 

4096 X 4096 

0.9 

2 

9% 

474% 

B-2 

0.09 

4096 X 4096 

0.9 

8 

1.9% 

100 % 

B-3 

0.09 

4096 X 4096 

0.9 

32 

0.42% 

22 % 

B-4 

0.09 

4096 X 4096 

0.9 

128 

0.05% 

2 .6% 

B-5 

0.09 

2048 X 2048 

0.9 

8 

4.2% 

212 % 

B-6 

0.09 

2048 X 2048 

0.9 

32 

0.9% 

45% 

B-7 

0.09 

2048 X 2048 

0.9 

128 

0.24% 

12 % 

B-8 

0.09 

2048 X 2048 

0.7 

8 

1 .6% 

80% 

B-9 

0.09 

2048 X 2048 

0.7 

32 

0.37% 

19% 

B-10 

0.09 

2048 X 2048 

0.7 

128 

0 .1% 

5% 

B-11 

0.09 

2048 X 2048 

0.5 

8 

0 .6% 

30% 

B-12 

0.09 

2048 X 2048 

0.5 

32 

0.13% 

7% 

B-13 

0.09 

2048 X 2048 

0.5 

128 

0.04% 

2 % 

G-1 

0.01 

4096 X 4096 

0.9 

2 

8.5% 

3000% 

G-2 

0.01 

4096 X 4096 

0.9 

8 

1.7% 

595% 

G-3 

0.01 

4096 X 4096 

0.9 

32 

0.36% 

126% 

G-4 

0.01 

4096 X 4096 

0.9 

128 

0.09% 

32% 

G-5 

0.01 

4096 X 4096 

0.9 

512 

0.03% 

10 % 

G-6 

0.01 

2048 X 2048 

0.9 

32 

0.75% 

265% 

G-7 

0.01 

2048 X 2048 

0.9 

128 

0.19% 

67% 

G-8 

0.01 

2048 X 2048 

0.9 

512 

0.08% 

29% 

G-9 

0.01 

2048 X 2048 

0.7 

32 

0.28% 

102 % 

G-10 

0.01 

2048 X 2048 

0.7 

128 

0.08% 

28% 

G-11 

0.01 

2048 X 2048 

0.7 

512 

0.044% 

15% 

G-12 

0.01 

2048 X 2048 

0.5 

32 

0 .11% 

39% 

G-13 

0.01 

2048 X 2048 

0.5 

128 

0.035% 

12 % 

G-14 

0.01 

2048 X 2048 

0.5 

512 

0.014% 

5% 


TABLE 2 

List of simulation runs used to test numerical convergence. All the runs are for a = 25 and Lx x = GOOdi x 388di and 

WERE PERFORMED OVER A DURATION OJpet = 3000. NOTE kT^ jTTleC? IS THE INITIAL PLASMA TEMPERATURE NORMALIZED BY REST ENERGY 

rriec^ . Time step is represented by the dimensionless Courant number Cr = cAt/Ar, where 
Ar = Ax Ay Az/(Ax Ay + AyAz + AxAz). NPC represents the number of particle pairs per cell. Eerr ! Etotai represents the 

RATIO BETWEEN CHANGE OF TOTAL ENERGY COMPARE TO THE INITIAL TOTAL ENERGY. Eerr!EkO REPRESENTS THE RATIO BETWEEN CHANGE 

OF TOTAL ENERGY COMPARE TO THE INITIAL PLASMA KINETIC ENERGY. 
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Fig. 12 .— Several cases with kTo = 0.36meC^ with grid number 
2048 X 2048, Courant number Cr = 0.7 but different numbers of 
particles per cell from 8 to 512. 
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Fig. 13.— Several cases for kTo = O.OlmeC^ with grid number 
4096 X 4096 and Cr = 0.9 but different numbers of particles per 
cell from 2 to 512. 
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